/*
 * =====================================================================================
 *
 *       Filename:  kmeans.h
 *
 *    Description:  K-means algorithm
 *
 *        Version:  1.0
 *        Created:  2009年06月15日 09时56分47秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */

#ifndef KMEANS_H
#define KMEANS_H
#include "ncvector.h"
#include "ncmatrix.h"
#include "mathutils.h"
/**
 *
 */
enum Distance {Euclidean, Cityblock, Cosine};

/**
 * \brief K-means algorithm
 *
 */
inline NCvector<int> kmeans(const NCmatrix<double> &data, int KK, NCmatrix<double> &means)
{
	int i,j,k,kmin,nchg=1,row = data.row(), column = data.column();
	NCvector<int> count(KK,0);
	NCvector<int> assign(row,0);
	means.resize(KK,column);
	NCvector<int> sel(KK);
	double dmin, d;


	//random select the K sample as the center.
	for( k=0; k<KK; k++)
	{
		for( j=0; j<column; j++ )
		{
			means[k][j] = data[(row-1)/(KK-1)*k][j];
		}
	}
//	std::cout<<"here"<<std::endl;
	//e-step
	while(nchg>0)
	{
		nchg = 0;
		for( k=0; k<KK; k++ )
			count[k] = 0;
		for( i=0; i<row; i++ )
		{
			dmin = 9e99;
			for( k=0; k<KK; k++ )
			{
				for(d=0.,j=0;j<column;j++)
				{
					d += SQE(data[i][j]-means[k][j]);
				}
				if( d<dmin )
				{
					dmin = d;
					kmin = k;
				}
			}
			if( kmin != assign[i])
			{
				nchg++;
			}
			assign[i] = kmin;
			count[kmin]++;

		}

		//m-step
		for ( k=0; k<KK; k++ )
		{
			for ( j=0; j<column; j++ )
				means[k][j] =0.;
		}
		for ( i=0; i<row; i++ )
		{
			for( j=0; j<column; j++ )
			{
				means[assign[i]][j] += data[i][j];
			}
		}

		for ( k=0; k<KK; k++ )
		{
			for( j=0; j<column; j++ )
			{
				means[k][j] /= count[k];
			}
		}

	}
	return assign;

}
#endif
